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NONLINEAR PHASE UNWINDING OF FUNCTIONS 


RONALD R. COIFMAN AND STEFAN STEINERBERGER 


Abstract. We study a natural nonlinear analogue of Fourier series. Iterative Blaschke factor¬ 
ization allows one to formally write any holomorphic function F’ as a series which successively 
unravels or unwinds the oscillation of the function 

F = aiBi + a2BiB2 + azBiB2B^ + ... 

where G C and Bi is a Blaschke product. Numerical experiments point towards rapid con¬ 
vergence of the formal series but the actual mechanism by which this is happening has yet to 
be explained. We derive a family of inequalities and use them to prove convergence for a large 
number of function spaces: for example, we have convergence in for functions in the Dirich- 
let space T>. Furthermore, we present a numerically efficient way to expand a function without 
explicit calculations of the Blaschke zeroes going back to Guido and Mary Weiss. 


1. Introduction 


1.1. Blaschke factorization. This paper studies a natural nonlinear way for unraveling the 
oscillation of a function F : C —)■ C that is holomorphic in a neighborhood of the unit disk. Our 
starting point is a fundamental theorem in complex analysis (Blaschke factorization) stating that 
any such function can be decomposed as 

F = BG, 


where F is a Blaschke product, that is a function of the form 


iei 


aj z- aj 
|ai| l-Wiz' 


where m G Nq and oi, 02 , • • • G D are zeroes inside the unit disk D and G has no roots in D. For 
1^1 = 1 we have |i?(z)| = 1 which motivates the analogy 


B ~ frequency and G ~ amplitude 


for the function restricted to the boundary. However, the function G need not be constant: it can 
be any function that never vanishes inside the unit disk. If F has roots inside the unit disk, then 
the Blaschke factorization F = B ■ G is going to be nontrivial (meaning B ^ 1 and G ^ F). G 
should be ’simpler’ than F because the winding number around the origin decreases and we will 
quantify this in many different ways. 


1.2. A formal series. There is a natural way of iterating Blaschke factorization that is inspired 
by the power series expansion of a holomorphic function in 0. Since G has no zeroes inside D, its 
Blaschke factorization is the trivial one G = 1 • G, however, the function G{z) — G(0) certainly 
has at least one root inside the unit disk D and will therefore yield some nontrivial Blaschke 
factorization G{z) — G(0) = BiGi. Altogether, this allows us to write 

F = BG 

= H.(G(0) + (G(z)-G(0))) 

= H.(G(0) + HiGi) 

= G{0)B + BBiGi. 
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At least formally, an iterative application gives rise to what we call the unwinding series 
F = Oj\B\ -\- CI2B1I32 (23-61-62-^3 A ^461626364 + ... 

This formal expansion first appeare(i in the PhD thesis of Michel Nahon [T 3 ]. Given a general 
function F it is not numerically feasible to actually compute the roots of the function; a crucial 
insight in m is that this is not necessary - one can numerically obtain the Blaschke product in a 
stable way by using a method that was first mentioned in a paper of Guido and Mary Weiss [ 25 ] 
(see also | 3 ]) and has been investigated with respect to stability by Nahon [T 3 ] and Letelier and 
Saito m- Numerical investigation m indicates that the formal series 

F = aiBi + 026162 + 03616263 + 0461626364 + ... 
will converge to the actual function and, generically, this seems to happen at an exponential rate. 

1.3. An example. The following example/picture is taken from the PhD thesis of Michel Nahon 
m- Let us consider the Blaschke factorization of a function given by the projection of a modulated 
Gaussian on the boundary onto holomorphic functions 

6(e*®) = 6+ . 

Fig. 1 shows the curves t —)■ 6(e*®), t —>■ 6(e*®) and t —>■ G(e®®) in the complex plane. A lot of 
the oscillation (and almost the entire phase) is transported from F to B leaving G significantly 
simpler than F. It also serves as a good example of the heuristic 

6 ^ frequency and G ~ amplitude. 

Figure 1 shows the real and imaginary part of the original signal, its shape when interpreted as a 
curve 6 : T —>■ C and the same information for 6 and G: 6 captures most of the oscillation. 


F on the real axis 


F in the complex domain 


Phase of F 



G on the real axis 



B on the real axis 







G (IGI > 0) in the complex domain 


B (IBI = 1) in the complex domain 


Phase of G 



Phase of B 



Figure 1 . A picture taken from Michel Nahon’s thesis [T 3 |: the behavior of the 
Blaschke decomposition of F on 9 D. 
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1 . 4 . Related work. Blaschke products have long been used in the signal analysis - often under 
the name Malmquist-Takenaka system. The crucial underlying fact is that for any two Blaschke 
products Bi,B2 the two functions Bi and B1B2 are orthogonal on L^( 5 D) since 

(Ri,BiB2)l2od) = / Ri(e*‘)Ri(e*‘)R2(e*‘)dt = / R2(e**)dt = 0. 

Jo Jo 

This allows naturally to build orthogonal functions via Bi, B1B2, BiB2B^, ... and contains the 
classical Fourier system l,z,z ^,... as a special case. We refer to papers of Eisner and Pap [ 5 ], 
Feichtinger and Pap [ 7 ], Pap m and Picinbino |14j for some examples. The unwinding series 
is first studied in the PhD thesis of Michel Nahon m- Subsequently, a method for numerical 
stabilization in the case of |F(e®®)| becoming small has been investigated by Letelier and Saito 
m- The unwinding series has been used by Healy [nmn] in the study of the Doppler effect. 
Of particular importance is a paper of Tao Qian |18j in which he proves the convergence of the 
unwinding series for F G 'H^(T). This paper was brought to our attention after this paper had 
been completed and we summarize his argument below. Closely related is also another approach 
developed by Qian, Ho, Leong and Wang [19] (and elaborated in further papers by Qian and 
collaborators (MllIilEolEg), which they call adaptive Fourier transform. The main idea is to use 
Blaschke products as a library and proceed by a projection pursuit approach, where at each step 
one projects onto the element in the library yielding the largest inner product with the function: 

fn+l = fn- {fn, Rn) Rn 

where is chosen among all Blaschke products with n zeroes as the one yielding the largest inner 
product. Since, in particular, the functions z" are elements of that library, this approach may 
be understood as a generalization of Fourier series - among their results is also an independent 
rediscovery of the Guido and Mary Weiss algorithm [T^ and of the unwinding series [ 511 . There 
are also similarities in spirit with recent work of Mallat [ 12 j . Mallat’s scattering transform is a 
translation-invariant operator, which is Lipschitz-continuous w.r.t. to C^— diffeomorphisms of 
the underlying space. The construction is based on an iterative application of wavelet transforms 
followed by restriction to the modulus. Our iterative application 

Gn{z) = G„( 0 ) -l- (G„(z) — G„( 0 )) = G„( 0 ) -I- R„+i(z)G„+i(z) 

uses the modulus of the corresponding functions while the coefficients are given as the mean. 
This yields a comparable level of stability: at least the leading coefficient is stable under both 
perturbations of the function and reparametrization of the torus. 

1.5. Notation and Outline. This paper deals with holomorphic ’signals’ given as functions 
/ : T —C by regarding them as the restriction of a holomorphic function F : C —>■ C on the 
boundary of the unit disk dD. We will therefore use both T and dD depending on which aspect 
should be emphasized. We will work with both Sobolev spaces FI and Hardy spaces "H (the 
Hilbert transform, which appears only briefly, will be T-Li). T) denotes the Dirichlet space, P+ the 
holomorphic projection. §2 states the results, §3 gives background material and discusses some 
possible applications. The proofs are given in § 4 . 

2 . Statement of results 

2 . 1 . Setup. Given a function F : C —>■ C, we define Gi as the outer part in the Blaschke factor¬ 
ization of F 

F = BiGi 

and then, iteratively, G„+i as the outer part in the Blaschke factorization 

Gn{z) — G„( 0 ) = Bn+i{z)Gn+l{z). 

We are interested in ensuring that ||G„||j(: —0 in some suitable space X and our main statements 
will be formulated that way. We emphasize that the formal series is, from the point of view 
of complex analysis, the canonical nonlinear extension of the Fourier series which arises from an 
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iterative application of G„(2;) — G„( 0 ) = z-Gn+i{z). The Blaschke series, in contrast to the Fourier 
series, proceeds by factoring out all zeroes inside D - this gives a rise to a much larger library of 
functions and makes it seem intuitive that one should not only expect convergence but also faster 
convergence than for the Fourier series. At the same time, the iteration 

Gniz) — G„( 0 ) = Bn+l(z)Gn+l{z) 

seems to define a very natural dynamical system on holomorphic function that could be of interest 
in its own right. 


2.2. Algorithm and roots. We start (assuming for simplicity that there are no roots on the 
boundary of the unit circle) with two basic observations. Recall that a general Blaschke factor 
has the form 




ai z — ai 
loil 1 - alz 


Suppose F : C —)■ C is holomorphic, F{0) ^ 0 and F has the set of roots R = {ri,r2,...}. If 
F = B ■ G, then the roots of G are simply 

{ Ti whenever \ri\ > 1 
1 /ri whenever 0 < \ri\ < 1. 

Geometrically, this means that the roots of F outside the unit circle stay unchanged while roots 
inside the unit circle are inverted across the unit circle. We emphasize that in every step of the 
algorithm consists of studying not G„(z) but G„(z) — G„( 0 ), which will have a very different set 
of roots. However, one immediate easy consequence is the following. 


Proposition. Let F : C —)■ C 6e given by a polynomial of degree n. Then the formal series 
converges and is exact after n steps. 

Proof. The algorithm is closed in the set of polynomials. Furthermore, since Gk{z) — Gfc(O) has 
at least one root in 0, the degree of the polynomial decreases by at least 1 in every step. □ 


This argument is more algebraic than analytic and comes with the obvious limitation that it does 
not give any convergence speed (analogously, it is not surprising that a trigonometric polynomial 
can be written as finite Fourier series). An illustrative example is given by 

F{z) = {z-{l-e)f-. 

The Blaschke factorization F = B ■ G is easy to write down and 

G(z) = (!-(!- £)z)^” = (1 - £)^” (z - . 

By making s sufficiently small, the functions F and G can be made as close to each other in any 
reasonable function space as we wish. These sort of examples immediately imply that it is not 
possible to construct a reasonable norm X with ||G„_|_i||jf < (1 — ^)||G„||x for some universal 
(5 > 0 . Exponential convergence, which is observed in practice, will therefore either not always 
be the case or be the consequence of an underlying phenomenon ensuring that iterative Blaschke 
factorization cannot always stay close to set of functions behaving like these polynomials. 


2.3. Regularity assumptions. Blaschke factorization only guarantees a splitting into an inner 
and an outer function, where the inner function itself is given by multiplying a Blaschke product 
and a singular inner function. It is not clear at this point how one would work with a singular inner 
function and we are restricting the further scope of the paper to functions that are holomorphic 
on a domain I? D D that contains an entire disk with radius I + e, where e > 0 can be arbitrarily 
small. This implies that the Blaschke factorization really factors into a Blaschke product and an 
outer product; moreover, any (nonzero) function that is holomorphic in a neighborhood of the 
unit disk has at most finitely many roots inside the unit disk, which guarantees that all Blaschke 
products are finite. This is not a serious restriction for applications as most signals of interest can 
be approximated by a trigonometric polynomial - it would be desirable to have a more complete 
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theory from a mathematical perspective, however, at this point even the dynamics of iterative 
Blaschke factorization on polynomials, though convergent, is far from being understood. 


2.4. A general contraction property. This section presents our main convergence result. We 
first state the result in the most general form and comment on special cases of particular interest 
further below. We start by introducing two norms on the Hardy space on the unit circle H^(D). 
Let 0 = 7o < 7i < ... be an arbitrary monotonically increasing sequence of real numbers and let 
X be the subspace of ■H^(T) for which 

2 



2 


E SD 

= 

E 

n>0 

X 

n>0 




< CXD. 


X 


n>0 


We define a second norm Y (semi-norm whenever 7 is not strictly increasing) 

2 



2 


E an 

= 

Y 

n>0 

Y 

n>0 




n>0 


Our main statement is that the Blaschke factorization acts nicely on these spaces. The first 
part of our statement is known (being ascribed to Digital Signal Processing in | 17 j l and can be 
equivalently phrased as follows: given a Blasche decomposition F = B ■ G and assuming both 
functions are expanded into a Fourier series 


n —0 n—0 

then, for every A G N 

00 00 

E \9n\^ < E 

n>N n>N 

Phrased differently, inner outer factorization shifts the energy to lower frequencies in a strictly 
monotonous way. Our main tool will be a refinement of that inequality. 


Theorem 1 (Main result). If F : D ^ C is holomorphic on some neighborhood of the unit disk 
and has a Blaschke factorization F = B ■ G, then 

iiG(e*')iiA < menwx- 

Moreover, if F{a) = 0 for some a G D, we even have 

||G(e*-)lll<l|J"(e")||i-(l-|an 

The most important implication is convergence of the unwinding series in the space Y if the 
initial data lies in X. The argument is straight-forward: the construction of the unwinding series 
proceeds by setting 

Bn+i{z)Gn+i{z) = G„{z) — G„( 0 ) 

and thus, by construction, the functions always have a root in a = 0 . Furthermore, adding and 
subtracting constants has no impact on because 70 = 0 and therefore 

||G„+i(e*')llr < WGnien - G„(0)||i - ||G„+i(e*')||^ 

= ||G„(e*')ll|-||G„+i(e*')||i 

Summing on both sides yields a telescoping series and thus 

00 

El|Gn(e*')llr<line*')llA, 

n=2 

which implies that ||G„(e®')||y —>■ 0 . After n steps, we have the equation 

F = OiBi + a2BiB2 ■ -f • • • B^-i + B1B2 ■ • ■ H„_i(G„ — G„( 0 )) 


G(e-) 


1 — az 
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and exploiting that \Bi \ = 1, we have that 

11/ — + • • • + Un-lBi ■ ■ ■ Bn-l)\\L'^{dB) = ||G„ — Gn (0) || ^2 (au). 

This motivates putting special emphasis on the space X arising from 7n = n for which Y = L^. 
This space is also known as the Dirichlet space T> and has special geometric significance and 
structure; for algebraic reasons we can get an even sharper inequality in that case (see below). 
Another natural (from a geometric perspective) space is given by 7 „ = where X = H^, 
Y = T> and Theorem 1 can be alternatively proven using Green’s formula (see below). All Sobolev 
spaces iJ® with s > 0 are also special cases: the statement implies that for A(e* ) G i?® with 
s > 0, we have convergence in iJ®“2. All these results have a completely analogous version on 
the upper half-space C+ with Blaschke-type products being defined on the real line M; even the 
proofs translate almost verbatim (see below). 

2.5. A slight generalization. The unwinding series can be phrased slightly more generally than 
we have done up to now: indeed, at the n—th step, we could actually pick an arbitrary q;„ € D 
and proceed via 

Bn+i(z) ■ G„+i(z) = Gn{z) — G„(a„). 

Clearly, Gn{z) — G„(q;„) is guaranteed to have at least one root in the unit disk because it has 
one in q;„ G D. Theorem 1 was formulated in a completely general way (for a general root a) and 
applies to this more general case as well (at the cost of introducing a factor 1 — |q;„P). In choosing 
a ‘good’ value for a„, one naturally encounters the quantity 

argmax(l — \z\^)\G{z)\. 

zeB 

We set an = 0, which - in practice - does not seem to make a big difference because the factor 
(1 — l^p) ensures that the maximum cannot be assumed on the boundary. While maximizing 
the quantity can lead to better results, we have observed that a„ = 0 seems to always be doing 
fairly well in practice. We will assume a„ = 0 throughout the rest of the paper but emphasize 
that the algorithm is slightly more general. One instance where this could be useful is whenever 
Gn{z) — G„(0) has a root on (9D: in order not to lose information on the phase, it is desirable for 
the performance of the Guido & Mary Weiss algorithm that |G„(z) — G„(0)| > 0 for all \z\ = 1. 
Whenever this is not the case, one could use G„(z) — G„(a) for a value a G D close to the origin 
such that this function has no roots on the boundary. 


2.6. A special case. Let us now explore the special cases with obvious geometric significance in 
greater detail. We identify functions F that are holomorphic in a neighborhood of the unit disk 
with maps 7 : T —>■ via 

7F(t) := F(e**). 

This is motivated by the fact that in the algorithm we obtain G„+i not from Gn{z) but from 
Gn{z) — Gn(0) and it is therefore natural to study translation-invariant (geometric) quantities 
depending on G„. Let us consider a particular example F(z) = (z + 0.3 + i/3)(z — 0.2)(z — 1.5 — i/2) 
(taken essentially at random) and the Blaschke factorization F(z) = B(z)G(z) (see Fig. 2 and 
Fig. 3). Since G{z) has no roots in D, the argument principle implies that G(z) does not wind 
around 0. Note that furthermore 

|F(e*‘)| = |G(e*‘)| 

for all t G R. As suggested by the picture (and many others like it), one would expect that the 
length of the curve jc is, at least generically, smaller than that of but we have been unable 
to prove that; instead, we were able to obtain that result for the natural version of length, 
sometimes called the energy of a curve 

p2-K 

energy(7) = / \-i' 

do 

By Holder’s inequality, we have that 


length(7)^ = 


p27r 


<27T 





2tt ■ energy(7). 
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Figure 2. F(e**) and G(e**) (dashed) for a cubic polynomial F. 


Therefore, in particular, if the energy of a curve tends to 0, then so will the length. Algebraic 
simplifications allow us to quantify the decrease of the norm of the boundary function in 
terms of its norm weighted against the Poisson kernel of the roots: the argument is not as 
sharp as the one formulated for the Dirichlet space further below but is very elementary (using 
Green’s theorem and geometric considerations). 

Theorem 2. Let F : D ^ C be holomorphic in some neighborhood of the unit disk. Then, if 
{aj : j G J} are the roots of F in D and F = B ■ G, 

p2-K p27z p27r _ -1 _ I |2 

/ |G'(e*«)pdt< / \F'{F 0 )\^dt- |G(e^*)PE | n~_ 

Jo Jo Jo Jgj ^jl 

Exploiting an additional geometric argument based on random projections and the uncertainty 
principle, we were able to obtain the following estimate, which controls the error in L°°(T). 


Corollary 1. Suppose F :£>—)■ C converges on some neighborhood of the unit disk. Then the 
formal series converges in L°°. Moreover, 

{n e N : ||G„( 2 ) - G„(0)|U=o(aD) > e}| < \FfF^)\^dt 

2.7. Winding numbers and the Dirichlet space. This section is entirely motivated by geo¬ 
metric considerations: we will discuss properties of closed curves in C given by 7 f(^) = F(e**). 
The winding number around a point z with respect to a curve 7 is defined as 

wind-y(zo) := A—: f - whenever zq ^ 7- 

2717 J.y z- Zo 

Examples strongly suggest that ’the average weighted winding number’ 



wind.Y(z)(iz 


should decrease. 


This quantity can be regarded as weighted area, which is the area enclosed by the curve weighted 
with the winding number. It arises naturally when one applies Green’s formula to compute the 
area surrounded by a simple, closed curve 7 : [ 0 , 27r] —>■ oriented counter-clockwise and written 
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as -fit) = {x{t),y{t)) via 


{x{t)y{t) - x{t)y{t))dt. 


Applying the very same formula in the case of a non-simple closed curve naturally gives rise to 


1 


{x{t)y{t) — x{t)y{t))dt = / wind-y(z)dz. 


This interpretation of the area formula dates back at least to a 1936 paper of Rado 
holomorphic, then we have 


]. If F is 


f windjp (z)dz = ( \F'{z)\'^dz. 

Jc Jn 



Figure 3. F(e**) (blue) given by the cubic polynomial. G'(e®‘) (dashed) has the 
same maximum winding but over a smaller area. 

Writing that representation in Fourier space gives the so-called area theorem stating that if 


f{z) = oo -I- aiz + a2Z^ -|-.... 

The Dirichlet space 

V = 


p OO 

then / wind.yj^ ( 2 )(i 2 : = tt n|a„ 

Jc 


/ : D —>• C| / holomorphic and / \f'{z)\^dz < oo i . 

was first introduced by Beurling and Deny [BE] .When equipped with the inner product 

{f,9)v = (/> 5)«2 + - / f'{z)g'{z)dz, 
jd 

it becomes a Hilbert space. A monotonicity statement for Blaschke decomposition in that space 
is well-known and follows at once from Carleson’s formula [5| (see also jS] Theorem 4.1.3]). 

Corollary 2 (Special case of Carleson’s formula). Assume F S 'H°°(D) with roots {aj : j € J} in 
D and has the Blaschke factorization F = B ■ G, then 


1 


\F'{z)\^dz= / \G'{z)\^dz+- 


/OD 




jCJ 


z - a, 


2 ■ 
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This result is better than Theorem 1 (which only gives the constant 1 instead of the sum over 
the Poisson kernel indexed by the roots) but follows from the same argument that we use to 
prove Theorem 1. This is due to some algebraic simplihcation that seems to only occur for 
X = T>,Y = and has to do with the fact that for a € D 

|z — q;| = |1 — az\ whenever \z\ = 1. 


2.8. A curious stability property. When doing Blaschke factorization F = B ■ G numerically, 
we will introduce some roundoff errors; even though we never actually compute the roots of the 
functions, this roundoff error can be imagined as perturbing the roots a little bit. We have the 
following curious and purely algebraic pointwise stability statement. 


Theorem 3. Suppose Fi,F 2 : C —>■ C are polynomials having the same roots outside o/D and the 
same number of roots inside D. Then the Blaschke factorizations 

Fi = BiGi and F 2 = B 2 G 2 , 


satisfy 


\Gi{z) — G 2 {z)\ = \Fi{z) — F 2 {z)\ whenever | 2 ;| = 1. 


This stability property was discovered by accident and seems quite curious. It is not clear to us 
whether there might be even more general statements of a similar type. 


2.9. An unwinding series on M. The inner-outer factorization was the crucial ingredient to 
our entire approach. A similar factorization can be achieved on the upper half-space. The role of 
Blaschke products is now played by functions indexed by Ai,..., A„ € C+ of the form 




which satisfies \B{z)\ = 1 on M. 


We will consider norms ll-l 


X > II iiy 


on the space 


T+(K) = {/ € : supp(7) C [0,oo)|. 


Let Ip : [0, 00 ] —>■ [0,oo] be a monotonically increasing, differentiable function with V'CO) = 0 and 


l/llx := / as well as 


11 - 




Theorem 4. If F has roots Ai,..., A„ G C+, then 




z - \k 


Z-\k 




X 


For the removal of a single root F(X) = 0, we have the stronger estimate 


F 


z — X 


<||F|||-(29(A)) 


X 


2 — AI 

Moreover, in the Dirichlet space ipif) = we even have 

|2 


X — X 


z — X 


Y 




z-Xk 


z- Xk 


X 


<\\Frx- [ \F{xrj2r^^dx, 


where the sum ranges over all roots of F on C+. 


3. Computation and application 

In this section we provide a collection of known facts, additional background material, a way of 
computing the Blaschke factorization without ever having to compute the roots (dating back to a 
1962 paper of Guido and Mary Weiss) and some sample applications. 
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3.1. Analytic signals. A classical way of using complex analysis when faced with a periodic, 
real signal u{t) : [0, 2t:\ —>■ K is to associate a natural imaginary part to the function. Already in 
1946 Gabor [5] argued that 

it has long been recognized that operations with the complex exponential 6-^“* [...] 
have distinct advantages over operations with sine or cosine functions. 


and proposed to analyze the signal 


f = u + iHiU instead, 


where Tti is the Hilbert transform. Vakman |24j proved that requiring certain natural assumptions 
on the complexification process, this is the canonical complexification. A convenient fact for actual 
computation is that if 

= ^ + ^ cos kO + bk sin 

k>l 


then 


{'Hiu){0) = ''^ak sin k9 — bk cos kO and {u + iHu){0) = ^ + (ufc — ibk)e^^^ 

k>l k>l 


3.2. The Guido and Mary Weiss algorithm. Let now f{9) be a complex signal (possibly 
obtained from a real signal using the process above). Assume additionally that f{9) ^ 0. Note 
that any such /(0) has only positive frequencies 

OO 

to which we may associate the function F : D —C 

OO 

= ^akZ^, 

k=0 

which, assuming sufficient regularity, has / as its boundary function. It is now our goal to construct 
the Blaschke decomposition oi F = B ■ G without computing the roots of the function. 

The algorithm proceeds as follows. 

(1) Compute the function g{9) = log \f{9)\. 

(2) Compute the analytic signal from g 

K^) = {9 + i'Hig){9). 

(3) Then we have the Blaschke factorization F = B ■ G, where 

= and B{9) = F{9)/G{9) 

on the unit circle. 

Clearly, the algorithm won’t work whenever there is a root on the boundary of the unit disk 
because log|/(0)| will be unbounded; also, whenever \f{9)\ becomes very small, the algorithm 
becomes unstable. Various ways for additional stabilization have been proposed: the stabilizing 
effect of adding a small constant has been investigated by Nahon m whereas Letelier and Saito 
m propose adding a small pure sinusoid. 


3.3. Removal of multiplicative noise. We return to the analogy 

B ^ frequency and G ~ amplitude. 

B is constructed from F by its roots inside the unit disk; conversely, given FiF') as boundary 
data, we can uniquely reconstruct the values of F inside D by convolving with the boundary data 
with the Poisson kernel. This compact integral operator enjoys a variety of smoothing properties; 
as a consequence it is stable against all sorts of perturbations (assuming they roughly preserve the 
local averages). A particular example given in (Fig. 4) consists of a function of the type 
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Figure 4. An example of a multiplicatively perturbed pure frequency cos (29) 
(dashed) on T and the resulting curve after two rounds of filtering. 


cos n6 ^ sinn0\ 

f{e) = cos (29) V an —^ + Pn—i=- , 

where a, j5 are instances of i.i.d. M(0, 1) random variables. We complexify the signal F and replace 
it by F{9)/\F{9)\. The outcome of two iterations of this process is shown in Fig. 4. A similar 
example can be found in work of Letelier & Saito m and Healy [siiin]. 


3.4. The instanteous phase problem. Given a complex signal, we may write it in polar coor¬ 
dinates as 


F(9) = |F(6»)|e*'^(®^. 


It is of interest in practice to understand how fast the frequency changes; naturally, if all quantities 
are well defined. 


F' = \F\'F‘f> + icl)'\F\F^ 


and thus 


F' 

~F 


1^1 




However, even assuming sufficient smoothness, the direct computation of the instaneneous fre¬ 
quency via this identity can be challenging and numerically unstable; various methods have been 
proposed (including one using Blaschke products due to Picinbono [II]). Blaschke products are 
well known to have the following particularly nice property: if the Blaschke product has finitely 
many roots, then 


imd 




and one has 



1 - IcifcP 




> m > 0. 


The unwinding series is therefore an approximation using strictly increasing frequencies, which 
greatly stabilizes numerical computation (see [ 13 ] for details). 


4. Proofs 

4.1. T. Qian’s Theorem. We start by giving a brief summary and proof of T. Qian’s theorem. 
This material is not new and can be found in |18j , however, that paper may not be easily accessible. 
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Theorem (T. Qian, IlSp. The unwinding series converges in for all 

OO 

F{9) = ^ with ^ |a„p < oo. 

n>0 n=0 

Proof from [18| . We first write the unwinding series in a slightly different way: since 0 is always 
a root in the iteration scheme, we may write it as 

F{z) = F{0) + jizBi + j 2 Z^BiB 2 + 732^^1 B2S3 + • • • + z^BiB 2 ... Bn{G - G(0)). 

We remark that any two of the Blaschke terms are orthogonal on L^(T): if £ < m, then 


e*‘) ) dt = 0 


j,2-!T / ^ \ / "* \ p2iT / m 

/ n n n 

•^0 \fe=i / \fe=i / -^0 \k=e+i 

because |i3fc(e®‘)| = 1 and the remaining term is holomorphic. We furthermore observe that 
the last term is orthogonal to all previous terms since the inner product simplifies by the same 
computation to 

•^0 \k=e+i / 

This immediately implies that 

+ ||e“*i?i(e“*)i?2(e“*) ■ • ■ B„(e*‘)(G(e*‘) - G(0))||^,(^^ . 

However, we can also guarantee that the remainder term is small by showing that it is orthogonal 
to all : 0 < fc < n — 1} since 

|^2^v / " \ 

7ne*‘” f n Bk{e^')] (^{e**) - G(0))e-*'=‘dt = 0. 

This implies 

OO 

||e™*Hi(e*‘)i?2 (e**)... H„(e*‘)(G(e**) - G(0))||^,(^) < \a,\\ 

k—n 

which then implies convergence as n —>■ cx). □ 


The proof shows that convergence will happen at least as quickly as Fourier series but potentially 
much faster since low-lying terms can already contain some part of the high-frequency contribu¬ 
tions. It would be interesting to quantifying how precisely this happens. 


4.2. Proof of Theorem We study the action of moving a single root from inside the unit disk 
D to the outside (inversion along the unit circle). Let |a| < 1 be the root; we compare 

f{z) = {z- a)F{z) and g{z) = (1 - az)F{z) 

on the boundary i9D. Expanding F\g^ into a Fourier series 

OO 

n=0 


we immediately get 


OO OO 

/(^)laD = “'^“0 + X! g(^)Ui[D = Qq + (q« - aan-i)z"- 

n—1 n—1 
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From the definition of || • ||x, we compute 

OO 

ll/Wlaoll^ =T'ol«n«oP + ^7n|an-i - aonP 

n—1 

OO 

— To I ^ I 1^0 I ^ ^ '~^n (I ^n— 11 1 “1“ |^| | O-n | ^ 

n=l 


and 

OO 

= Tolo-ol + ^ ^ 7n|^in ~ CK^n-ll 

n—1 

OO 

To I ^0 I ^ ^ Tn (I I ^(^n—l^n CXCtnClrL— 1 “1“ | ^ I |^n—11 ) ■ 

n—1 


We see that the mixed terms appear in both sums and cancel: subtraction yields 

OO 

ll/(^)|aDll^ “ llff(^)|a]D,llA = -7o(l - laHlaoP + (1 - \a\'^)^in{\an-i\^ - |a«T 

n—1 

OO 

= (1 - lap) 51 (^"+1 “ T'") l“"l^ 

n—0 

= (l-|ap)||F^|^„||^. 

This equation has a nice and definite form but we will only use it in one instance. Let us assume 
we are given F(Z) and a finite list of roots {ai, a 2 , ■ ■ ■, an} C D. We know, by construction, that 
at least one of the roots is 0 and we assume without loss of generality that a„ = 0. Then we can 
consider the sequence of functions 


F{z) = {z- ai){z - a2){z - as)... (z - a„_i)(z - a„)n(z) 
^i(^) = (1 - Wz)iz - a2){z - Us)... (z - a„_i)(z - an)H{z) 
F2 {z) = (1 - aTz)(l - a^z)(z - as)... (z - q;„_i)(z - an)H{z) 


Fn-i{z) = (1 - q;iz)(1 - q;22)(1 - a^z) ... (1 - a„_iz)(z - an)H{z) 
and we can conclude from the computation that 

\\m\oJ\x > \\Fi{z)\J\x >■■■> ||F„_i(z)|^^||x. 

Clearly, the outer function G in the Blaschke decomposition F = B ■ G will be given by 
G'(z) = (1 - alz)(l - a^z)(l - aifz)... (1 - aFriz)(l - ^z)H{z) 

In the final step, we use the fact that there is always one root satisfying q;„ = 0 and exploit the 
full strength of the argument to conclude that 


F„_i(z)|^„|||-||F„(z)^J|| = 




P[ (1 - akz)H{z) 


an 


k=l 


= l|G(z) 


Y 


More, generally, if there is no root in 0, then applying the same argument yields 


F„_i(z)|^„||| - ||F„(z) = (1 - |a„P) 


aB" 


P[ (1 - akz)H{z 


lao 


fc=i 


= (i-KP) 


Gjz) 

1 - oFz 


I an 


This concludes the argument. □ 


y 
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The last part of the argument highlights a fundamental difficulty: while there is an effective gain 
every time we move a root to the outside, it is not clear to us how the sum of these gains could 
be properly controlled (which is why we only take the last one). This we only managed to do in 
the case of the Dirichlet space, where an additional (algebraic) simplification takes place. 

4.3. Proof of Theorem]^ We study again the action of moving a single root to the outside by 
inversion along the unit circle. The computation resembles the computation in the more general 
case except that we are able to invoke Green’s formula at the end of the argument. 

Lemma 1. Let F be analytic in a neighborhood of the origin and a G C with |a| < 1. If 

f = {z — a)F and g = {1 — az)F 

then 

p‘2'k p‘2'k p2TT 

/ \g'{F<^)fdt< \f{F<^)\^dt-{l-\af) \F{F*)fdt 

Jo Jo Jo 

whenever all terms are defined and finite. 


Proof. Obviously 

f = F + {z- a)F' 

and thus 

|/'|2 = |f|2 + F{z - a)F' + F{z - a)F' + |z - a|VT- 

At the same time 

g' = -aF + (1 - az)F' 

|g'|2 = \a\^\F\^ - aF{l - az)F' - aF{l - az)F' + |1 - az\'^\F'\^. 

If \z\ = 1, then jz — a|2 = |1 — azp and since we only integrate over 9D, we get 

[ |/'|2_ = (i_|a|2) /■ \F\^+ f F{z - a)F' + F{z - a)F' + aF{l - az)F' + aFjl 

Jda Job Jdo 

This is already almost what we want, it remains to show that 


dz)F'. 


Ida 


F{z - a)F' + F{z - a)F' + aF{l - az)F' + aF{l - az)F' > 0. 


The expression can be rewritten as 


which is 


f FF'{{z — a) + a(l — az)) + FF'{[z — a) + a(l — az)), 
JdB 

(l-lap)/ FPz + FF'z, 


/9D 


Now we go back from the classical derivative f'{z) to the angular derivative along the boundary 
of the disk f{z). As before 

which can be rewritten as 

f'[z) = —izf(z) whenever \z\ = 1. 

Using this, we can rewrite the terms as 


FzF' = Fz{—i)zF = iFF whenever \z\ = 1 
F'~zF = —izzFF = —iFF whenever \z\ = 1. 


We need to show that 


FF-FF> 0. 


IdD 


If we write 


U(e**) = x{t) + iy{t) 











NONLINEAR PHASE UNWINDING OF FUNCTIONS 


15 


then 

i{FF - FF) = 2{x{t)y{t) - x{t)y{t)). 
The problem consists now of evaluating 


/ 2(x(t)y(t) - 

J dn 

This corresponds to integrating the vector field 

2{—y,x) over the curve 
Green’s theorem states that this implies 


x{t)y{t))dt. 


lit) = ixit),yit)). 


2{x{t)y{t) 


x{t)y{t))dt = 4^ > 0, 


JdB 

where A is the area of the domain enclosed by the curve 7 (weighted at each point with the winding 
number with respect to 7). □ 


li F = B ■ G has more than one root in D, Lemma 1 can be applied iteratively. 


Proof of Theorem 2. The previous language establishes a relationship between /, g and F. How¬ 
ever, only the modulus of F ever appears in the argument: exploiting that 

\z — ai \ = |1 — <xiz\ whenever \z\ = 1 

allows for a better using of the gain obtained from iterative application of the previous Lemma 
when inverting several roots along the unit circle. More precisely, consider again 

F(z) = {z- ai){z - a 2 )iz - as) ...{z- a„_i)(z - a„)iL(z) 

Fi(z) = (1 - aiz){z - a 2 )iz - as) ...{z- a„_i)(z - a„)H{z) 

F 2 iz) = (1 - aiz){l - a^z)(l - cdsz) ...{z- q;„_i)(z - an)H{z) 


Fn-i{z) = (1 - aiz){l - a 2 z)(l - asz) ... (1 - an-iz){z - an)H{z) 

The crucial new ingredient is that 

ll^(^)|aDll-£'’'(9») = = • • • = ll^n-l(^)|aiD,llL2(9D)- 

The very same reason allows for a more precise analysis of the effect removing one root has. Let 
again |q;| < 1 be the root; we compare 

f{z) = iz- a)F{z) and g{z) = (1 - az)F{z) 

on the boundary i9D. The same computation as before yields 

ll/W|aDlll - hiz)\oj,\\l = (1 - \c^\")\\fU\1^ = (1 - i«n 

In particular, all the arising expressions can be summed in closed form and the arising gain is 

mil - iiGiii = iGLUi-y; 

□ 

There is a difference of a factor 2 in the way we stated the Carleson’s formula and the proof 
above: this is due to the fact that we computed the effect on what turns out to be whereas 
the Dirichlet space in Carleson’s formula also has the norm (which stays preserved since 
|F| = |G|), hence the difference of a factor of 2. 
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4.4. Proof of Corollary We have 

P‘2'K p‘2.'k p 

/ / |G;{e*^)|2d0- / |G„(e*«)-G„(O)pd0, 

Jo Jo Jdn 

but there is no way of turning this into a quantitative decay estimate because the gain 


|G„(e*®)-G„(O)|^d0 


need not be proportional to the size of 



|G;(e*«)pd0. 


Put geometrically, G„(e*®) may wind around G„(0) very quickly while |G„(e*®) — G„(0)| could 
be quite small all the time. The crucial insight is as follows: if that were actually the case and 
||G„(e*®) — G„(0)||i2(T) is small, then one would certainly hope that |jG„(e*®) — G„(0)||iooc]r) is 
also small. Now we reverse the order of that argument: suppose that ||G„(e*®) — G„(0)||ioo(T) is 
not small. This means that |G„(e*®) — G„(0)| is big for some 9, which does not at all mean that 
the function is large in it could just be big in one place and very small everywhere else: this, 
however, would imply that the norm of the gradient is large and we know it cannot exceed 
that of the initial data. 


4.4.1. Sobolev embedding. The second ingredient of the argument may be formulated as follows: 
let h : T —)■ K be differentiable. If ||h'||j^2(T) is not very big and ||h||iDo(T) has a certain size, then 
II^IIl 2 (t) cannot be arbitrarily small (depending on the first two quantities): the only way to be 
big in L°°(T) but small in L^(T) is quick decay around the point where the supremum is assumed. 
The inequality is merely the classical embedding of the Sobolev space i?^(T) L°°(T). 


Lemma 2. Let h :T ■ 


be a differentiable function which changes sign. Then 


Il 2 > 


16 \\h' 


t 


Proof. Assume without loss of generality that h{0) = 0. Assume x to be such that \h{x)\ 
Using the Cauchy-Schwarz inequality, we get 

\\h\\l^(T) = \h{x)\^ \h{z)\^dz^ \h'{z)\^dz^ . 

Squaring both sides gives the result. 


□ 


4.4.2. Random projections. We apply the statement to a curve : 9D —>■ C, which is different 
object than a periodic function h : T —>■ K. The natural approach to reduce one to the other would 
be to fix a vector n € with unit length \v\ = 1 and consider the projection 

Ki) = (7(0.■ 


The next Lemma states that there exists a unit vector v such this reduction does not change the 
norm and norm by more than an absolute constant: 

P‘27Z p^ItZ i 

\lit)\^dt<6 \{n,^{t))\^dt. and |h(0)| = | (7(0),:/) | > ^|7(0)|. 

It is easy to see by Cauchy-Schwarz that h varies slower than 7 



i(y(o,^)i<i7'(oi- 


Therefore, after having established the existence of such a vector n and reparametrizing the curve 
in such a way that |/i(0)| > 2“^/^ maxj |7(t)|, we could deduce that 


\l{t)'^dt> / \h{t)\^dt> 


1 11/^11 


16 \\h' 


li= 


> 


1 maxt|7(t)0 

64 W{z)\Hz 


which is a quantitative version of our intuition described above: in order for the function to be big 
at some point, it cannot be too small on average. Let us now prove the statement. The argument 
says that it is sufficient to take that vector at random to have the desired property to be true on 
average (in particular, there exists at least one vector for which it is true). 
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Figure 5. Assuming arclength-parametrization, projecting the curve onto the 
dashed line preserves both norm and the norm up to absolute constants. 


Lemma 3. Let 7 : T —>■ be a periodic curve in the plane and assume 7(0) ^ (0, 0). Then there 

exists a unit vector lul = 1 with 


(7(0),u) I > ;;^l7(0)| 


as well as 


r2n 


r27T 


< 6 


(z/,7(t)) \^dt. 


Proof. The line from the origin to 7(0) defines a unique angle a. Let us now chose v randomly 
from a — 7r/4 and a + 7r/4. Any such vector satisfies the first condition and we will now compute 
the expectation of the norm for such a random vector. We first remark that for every fixed 
vector ^ S and every 0 < a < 27r a simple computation shows that 

/■“+5 r\ 

/ I ((coss,sins),£) / | ((coss,sins), (0,1)) pds 

> f (sins)^ ds 


= 14 


, TT — 2 


We now compute the expectation by using this and exchanging the order of integration 


2 

TT 




q-7 Jo 


f.27r /-a+f 


((coss,sins),7(t)) \^dtds = — 

77 , 


((coss,sins),7(t)) \^dsdt 


0 Ja-f 

271- 


2 TT- 2 


> 


TT 4 


mfdt. 


Since (27r)/(7r — 2) < 6 and since a random vector has that expectation, there exists at least one 
vector with that value. □ 

4.4.3. Proof of Corollary^ The monotonicity formula implies that 


^271- 


(.2-7r 


|(G„(e®') — G„(0))'fdt = / |G(j(e*‘)fdt is monotonically decreasing in n. 


Suppose that for some n and some z G 


We can identify G„{z) — G„(0) : 
Lemma so that 


|G„(z)-G„(0)| >£. 
A C with a curve 7 : T 
e 


and reparametrize it using our 


(7(0), u) I > 


42 
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and 


Jo Jo 

We can now apply our second Lemma to the function 

h(t) = ■ 

Since G„(z) — G„(0) has winding number at least 1, so has 7 and therefore h vanishes at least 
in two points, h and Gn{z) — G„(0) have comparable L°°(T)—norm (up to a factor of ^/2) and 
comparable L^(T)—norm (up to a factor of 6) and elementary geometric considerations (projections 
make vectors only smaller) show that the derivative of h satisfies 

\h'{t)\ < |7'(t)|. 

We can now apply Lemma 2 and conclude that 


|(G„(e*‘)-G„(0))|'dt^ 




lli= 


> 


11^1112 /y |i^'(e*«)|2d6»’ 


where the second inequality follows from the assumption that |G„(z) — G„(0)| > e and the fact 
that the L'^—noim of the derivative is decreasing. Now, let’s look at the next step in the algorithm, 
where we decompose 

Gn — Gn(0) = BGn+l- 

Our inequality tells us that the squared L'^—noru\ of the derivative decreases at least by a factor 
of (using \B\ = 1 for Blaschke products) 


n27T p2TT 

/ |G„+i(e*‘)pdt= / |G„(e**)-G„(0)pdt> 

Jo Jo 




This yields 


r2Tz 


^277 


^277 


|G;+i(e*«)pdt< / \G'^{e^0)\^dt- |G„(e“‘) - G„(0)|2dt 


< 


p2-K 

/ |GUe^^)pdt-c 
Jo 


/y \F'{e^^)\He 


for some universal constant c > 0. However, since all the involved quantities are nonnegative, this 
immediately implies that number of n G N for which 

|G„(z) — G„(0)| > e for some z G i9D 

is bounded from above by 


\F'{F^)\^d0^ / l^c- 






This concludes the argument. □ 


4.5. Proof of Corollary We recall the action of removing one root which entails comparing 
f{z) = {z-a)F(z) and g{z) = {I - az)F{z). 

As was shown in the proof of Theorem we have the identity 

ll/W|aDll|-|l5W|aDllx = (l-|«P)ll^laDll^- 

In the Dirichlet space X = T), we have 7„ = n and thus Y = L?. In the proof of Theorem 1, we 
used the monotonicity formula to remove all roots and applied the full strength of the inequality 
only once for a root that is in the origin. Here, the special algebraic structure of the space allows 
us to apply the inequality multiple times and sum all the contributions in closed form. The crucial 
ingredient that makes this possible is the algebraic identity on 9D for all |a|, |/3| < 1 

\z — a\\z — I5\ = |1 — az||z — /3| = |1 — az||l — f3z\ whenever \z\ = 1. 
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We will now illustrate the effect of applying the identity twice (to remove two roots a, /3 from D). 
The arising functions are 

fi{z) = {z-a){z-/3)F{z), f 2 iz) = {l-az){z-l3)F{z) and fsiz) = (1-az){l-j3z)F{z). 
Applying the identity twice yields 

ll/i Wlanll^ - ll/2Wlaell^ = IK" - 

\\f2{z)\J\l-\\hiz)\J\l = |i(l-az)^^|,„||i. 

Normally, we would be unable to sum up these two contributions, however, here the algebraic 
identity implies that 

and 

Ul-az)F\J\l, = \\{z-a)F\gJ\l.= 

This allows us to simplify 

ll(» - «F|*lli- + lO - 5*)F|a„lli. = 

The general case for more sums follows by the same reasoning. □ 


fi{z) 


z - . 


I9D 


L2 


fi{z) 


z — a 




12 


The argument can be easily summarized as saying that the algebraic structure oi X = V implies 
that Y = L^-, the additional algebraic ingredient is |i?(e*®)| = 1 which implies that the various 
||F IooIIt one gets from successive removal of roots can actually be summed up in closed form. 

4.6. Proof of Theorem [sl 


Proof. The statement is pointwise and invariant under multiplication with polynomials having all 
roots outside of D: it thus suffices to prove it for polynomials having all their roots inside of D. 
We write 


Obviously 


and thus 


Fi=Y[{z-ai) and F 2 = Y]_iz - I3i). 




i-1 


i?i = TT ^— — as well as ^2 = TT —= 

^ 111 - TWt: 11 1 _ 


Z- I3i 


i=l 


1 - aiZ 


i=l 


1 - &,.Z 


Gi = (1 — aiz) as well as G2 = (1 — Piz). 


i=l 


i=l 


An explicit computation yields that 


gi(.)-g2(z)=^.m ^ n(-^)- E n(-K.) 


k=0 \ AC{l,...,n} j^A 


jeB 

B|=fc 


while 


f,{z)-f2{z)=j2z'^\ E • 


k=0 \ AC{l,...,n} jl^A 

' |A|=r.-fe 


sell ,jeB 


Altogether, this implies that if 

n n 

Gi{z) - G 2 {z) = then Fi{z) - F 2 {z) = y^ a^z^. 




/c=l 
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It remains to show that both quantities have the same norm if |z| = 1 


\F,{z)-F2{z)\ = 


^C„_ 


kZ 






kZ 




%—k 




—n—k 

k z 


/c=l 


^ ^ ^n—kZ 


k^l ^ 

* = |Gi(z) - G2iz)\. 


□ 


4.7. Proof of Theorem [31 


Proof. We imitate the argument in the case of Fourier series and again study the effect of removing 
one root by comparing 


Note that 


and therefore 




Ia-II5IIa = 


After simple computation we arrive at 


- a)F{z) 

and 

9{z) = 

- 

and 

II 

^ d - 


^ . d 


- 

~ 




ipi^d^ 


t-F{0-aF{0 


t-F{0-aF{0 


2 

r — —1 

= i{a — a) 

F'F + FF' 


Writing P(^) = a(^) + z&(^) as real and imaginary parts, we see that 

F'F + FF' = 2{a'a + b'b). 


This implies that we can write 


d - 

2 


, F{0-aF{i) 

df, 

— 

, F[0-c.F{0 

df, 


and therefore with integration by parts 


^(a-a)^\Fm^ 

d^ 


WfWx - Mx = (*(« - = 25a \F{0?^'{0d^, 


which is clearly nonnegative because a S C+. It remains to study the special case of the Dirichlet 
space: if ^ we have with Plancherel that 


poo poo p 

/ |F(oivm=/ \f{om= m)\^dc 

Jo Jo Jr 

The key ingredient is again of an algebraic nature: the difference can be quantified in terms of a 
quantity whose behavior can be controlled while removing several roots one after the other. Let 
us illustrate this again with 


fi{z) = (z - a){z - P)F{z), f 2 {z) = {z-a){z-ld)F{z) and f^iz) = {z - a){z -/3)F{z). 


Applying the identity twice yields 

\\fl{z)\\l-\\f 2 {z)\\l = \\{z-m\h 
\\f2iz)\\l-\\Mz)rv = \\iz-a)F\\l. 

and we see once more that the sum of the gain can be controlled. This yields the statement. □ 
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